{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Gromov-Wasserstein Barycenter example\n",
    "\n",
    "\n",
    "This example is designed to show how to use the Gromov-Wasserstein distance\n",
    "computation in POT.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Erwan Vautier <erwan.vautier@gmail.com>\n",
    "#         Nicolas Courty <ncourty@irisa.fr>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "\n",
    "import scipy.ndimage as spi\n",
    "import matplotlib.pylab as pl\n",
    "from sklearn import manifold\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "import ot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Smacof MDS\n",
    "----------\n",
    "\n",
    "This function allows to find an embedding of points given a dissimilarity matrix\n",
    "that will be given by the output of the algorithm\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def smacof_mds(C, dim, max_iter=3000, eps=1e-9):\n",
    "    \"\"\"\n",
    "    Returns an interpolated point cloud following the dissimilarity matrix C\n",
    "    using SMACOF multidimensional scaling (MDS) in specific dimensionned\n",
    "    target space\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    C : ndarray, shape (ns, ns)\n",
    "        dissimilarity matrix\n",
    "    dim : int\n",
    "          dimension of the targeted space\n",
    "    max_iter :  int\n",
    "        Maximum number of iterations of the SMACOF algorithm for a single run\n",
    "    eps : float\n",
    "        relative tolerance w.r.t stress to declare converge\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    npos : ndarray, shape (R, dim)\n",
    "           Embedded coordinates of the interpolated point cloud (defined with\n",
    "           one isometry)\n",
    "    \"\"\"\n",
    "\n",
    "    rng = np.random.RandomState(seed=3)\n",
    "\n",
    "    mds = manifold.MDS(\n",
    "        dim,\n",
    "        max_iter=max_iter,\n",
    "        eps=1e-9,\n",
    "        dissimilarity='precomputed',\n",
    "        n_init=1)\n",
    "    pos = mds.fit(C).embedding_\n",
    "\n",
    "    nmds = manifold.MDS(\n",
    "        2,\n",
    "        max_iter=max_iter,\n",
    "        eps=1e-9,\n",
    "        dissimilarity=\"precomputed\",\n",
    "        random_state=rng,\n",
    "        n_init=1)\n",
    "    npos = nmds.fit_transform(C, init=pos)\n",
    "\n",
    "    return npos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data preparation\n",
    "----------------\n",
    "\n",
    "The four distributions are constructed from 4 simple images\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rflamary/.local/lib/python3.6/site-packages/ipykernel_launcher.py:6: DeprecationWarning: `imread` is deprecated!\n",
      "`imread` is deprecated in SciPy 1.0.0.\n",
      "Use ``matplotlib.pyplot.imread`` instead.\n",
      "  \n",
      "/home/rflamary/.local/lib/python3.6/site-packages/ipykernel_launcher.py:7: DeprecationWarning: `imread` is deprecated!\n",
      "`imread` is deprecated in SciPy 1.0.0.\n",
      "Use ``matplotlib.pyplot.imread`` instead.\n",
      "  import sys\n",
      "/home/rflamary/.local/lib/python3.6/site-packages/ipykernel_launcher.py:8: DeprecationWarning: `imread` is deprecated!\n",
      "`imread` is deprecated in SciPy 1.0.0.\n",
      "Use ``matplotlib.pyplot.imread`` instead.\n",
      "  \n",
      "/home/rflamary/.local/lib/python3.6/site-packages/ipykernel_launcher.py:9: DeprecationWarning: `imread` is deprecated!\n",
      "`imread` is deprecated in SciPy 1.0.0.\n",
      "Use ``matplotlib.pyplot.imread`` instead.\n",
      "  if __name__ == '__main__':\n"
     ]
    }
   ],
   "source": [
    "def im2mat(I):\n",
    "    \"\"\"Converts and image to matrix (one pixel per line)\"\"\"\n",
    "    return I.reshape((I.shape[0] * I.shape[1], I.shape[2]))\n",
    "\n",
    "\n",
    "square = spi.imread('../data/square.png').astype(np.float64)[:, :, 2] / 256\n",
    "cross = spi.imread('../data/cross.png').astype(np.float64)[:, :, 2] / 256\n",
    "triangle = spi.imread('../data/triangle.png').astype(np.float64)[:, :, 2] / 256\n",
    "star = spi.imread('../data/star.png').astype(np.float64)[:, :, 2] / 256\n",
    "\n",
    "shapes = [square, cross, triangle, star]\n",
    "\n",
    "S = 4\n",
    "xs = [[] for i in range(S)]\n",
    "\n",
    "\n",
    "for nb in range(4):\n",
    "    for i in range(8):\n",
    "        for j in range(8):\n",
    "            if shapes[nb][i, j] < 0.95:\n",
    "                xs[nb].append([j, 8 - i])\n",
    "\n",
    "xs = np.array([np.array(xs[0]), np.array(xs[1]),\n",
    "               np.array(xs[2]), np.array(xs[3])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Barycenter computation\n",
    "----------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ns = [len(xs[s]) for s in range(S)]\n",
    "n_samples = 30\n",
    "\n",
    "\"\"\"Compute all distances matrices for the four shapes\"\"\"\n",
    "Cs = [sp.spatial.distance.cdist(xs[s], xs[s]) for s in range(S)]\n",
    "Cs = [cs / cs.max() for cs in Cs]\n",
    "\n",
    "ps = [ot.unif(ns[s]) for s in range(S)]\n",
    "p = ot.unif(n_samples)\n",
    "\n",
    "\n",
    "lambdast = [[float(i) / 3, float(3 - i) / 3] for i in [1, 2]]\n",
    "\n",
    "Ct01 = [0 for i in range(2)]\n",
    "for i in range(2):\n",
    "    Ct01[i] = ot.gromov.gromov_barycenters(n_samples, [Cs[0], Cs[1]],\n",
    "                                           [ps[0], ps[1]\n",
    "                                            ], p, lambdast[i], 'square_loss',  # 5e-4,\n",
    "                                           max_iter=100, tol=1e-3)\n",
    "\n",
    "Ct02 = [0 for i in range(2)]\n",
    "for i in range(2):\n",
    "    Ct02[i] = ot.gromov.gromov_barycenters(n_samples, [Cs[0], Cs[2]],\n",
    "                                           [ps[0], ps[2]\n",
    "                                            ], p, lambdast[i], 'square_loss',  # 5e-4,\n",
    "                                           max_iter=100, tol=1e-3)\n",
    "\n",
    "Ct13 = [0 for i in range(2)]\n",
    "for i in range(2):\n",
    "    Ct13[i] = ot.gromov.gromov_barycenters(n_samples, [Cs[1], Cs[3]],\n",
    "                                           [ps[1], ps[3]\n",
    "                                            ], p, lambdast[i], 'square_loss',  # 5e-4,\n",
    "                                           max_iter=100, tol=1e-3)\n",
    "\n",
    "Ct23 = [0 for i in range(2)]\n",
    "for i in range(2):\n",
    "    Ct23[i] = ot.gromov.gromov_barycenters(n_samples, [Cs[2], Cs[3]],\n",
    "                                           [ps[2], ps[3]\n",
    "                                            ], p, lambdast[i], 'square_loss',  # 5e-4,\n",
    "                                           max_iter=100, tol=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualization\n",
    "-------------\n",
    "\n",
    "The PCA helps in getting consistency between the rotations\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fa5bf8b5cc0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAJDCAYAAABHZBNLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3V2MJdV5//vfM2Agzf8ohgFhbJge5ogYYyUnMS3LxFLkExMJcwF27ERYLQKRUcfYlq/HGulEsjSKk8vIVjgtFEFCywYjWR4nRMiv8kUEdiMZ83YwA2LGQ8Ae8F+WrIlw7HnORVUzu/fs2vW2qtaqqu9HWur9Ul2revYzaz211qra5u4CAABAHHtiHwAAAMCUkYwBAABERDIGAAAQEckYAABARCRjAAAAEZGMAQAARBQkGTOzfzazn5vZUwXvm5n9o5kdNbMfm9l7QtSL8SCGEAJxhLaIIcQQamTsXkk3Lnn/Q5KuzsuGpH8KVC/G414RQ2jvXhFHaOdeEUPoWZBkzN2/L+kXSza5RdK/eOZRSW81s8tD1I1xIIYQAnGEtoghxNDXmrF3SPrpzPMT+WtAVcQQQiCO0BYxhODOjX0As8xsQ9mwry688MLrrrnmmshHhK49/vjjr7n7pSH3SRxNSxcxJBFHU0NbhLbaxFBfydjLkq6ceX5F/tou7r4paVOS1tbWfHt7u5+jQzRmdqzippViSCKOpqZGDEnEEQrQFqGtmm3RLn1NUx6R9Ff5VSjvk/RLd3+lp7oxDsQQQiCO0BYxhOCCjIyZ2ZclfUDSJWZ2QtLfSnqLJLn73ZIelnSTpKOSTkn66xD1YjyIIYRAHKEtYggxBEnG3P3jJe+7pE+HqAvjRAwhBOIIbRFDiIE78AMAAEREMgYAABARyRgAAEBEJGMAAAARkYwBAABERDIGAAAQEckYAABARCRjAAAAEZGMAQAAREQyBgAAEBHJGAAAQEQkYwAAABGRjAEAAEREMgYAABARyRgAAEBEJGMAAAARkYwBAABERDIGAAAQEckYAABARCRjAAAAEZGMAQAARBQkGTOzG83sOTM7amYHF7x/h5mdNLMf5eXOEPViXIgjtEUMIQTiCH1rnYyZ2TmSviTpQ5KulfRxM7t2waYPuPsf5uWetvViXIijMLa2pP37pT17sp9bW7GPqD/EEEIgjmqacqMTUIiRsfdKOuruL7r7ryV9RdItAfaLWeMPeOKopa0taWNDOnZMcs9+bmyMMVQKEUMIgTiqikYnmBDJ2Dsk/XTm+Yn8tXkfNbMfm9lDZnZlgHqnYxoBTxy1dOiQdOrU7tdOncpenwhiCCEQR1XR6ATT1wL+b0ja7+5/IOmbku5btJGZbZjZtpltnzx5sqdDGwACfgdxtMTx4/Ven6hKMSRNN45QCW2RRKMTUIhk7GVJs2cFV+SvvcndX3f3N/Kn90i6btGO3H3T3dfcfe3SSy8NcGgjMY2AJ45qWDRrvW/f4m2LXh+hYDGUbzv6OMJCtEVVl8XQ6AQTIhn7oaSrzewqMztP0q2SjsxuYGaXzzy9WdKzAeodn6L/ANMIeOKooqJZ65tuklZWdm+7siIdPhznOCMghhIwguWt046jOstiDh+efKMTjLu3LpJukvQTSS9IOpS/9nlJN+eP/07S05KekPRdSdeU7fO6667zSbn/fveVFfcs/LOyspK9vuy9gZO07cRRLauru0Nhp6yuZiGxuupudub52HUdQz6gOIr9+Q+5qaItyi1rYBaJHXQJmY2husWy30/P2tqab29vxz6M/uzfn52BzFtdlV56KTsrOXQom5rcty8781hf7/sogzOzx919rav9jzGO9uzJWsd5ZtLp0/0fT2xdx5A0jDjaGdCYXV66siJtbvbXVJQ1YymjLcrRwDTWJoa4A38qytaFra9nrdnp09nPESRiaKburPUIpo1QQQrX+UxjeevIlTUwNCidIBnr07IgpodFRXWWaUzjriiQ0kiEprG8deSWNTA0KJ0hGetLWRDTw6Ki9fVs6ml1NZs5WF0tnopKYbQE/YiZCO2cGx47lsXkLNZzD8yyBoYGpTMkY30pC2J6WNSwM2v9r/+aPb/ttsUDpCmMlqAfXV3YVjYIP3tuKGXnhzsJ2bJmDAkrWhZTpUFh1qYRkrG+VAni+f8A0uKgpoeFqg2QMm00HXXO56qqEmOLzg3dzyzaJxEbkSrryZi1aYRkLLRQ9wpbFtT0sFC1AVJuAzQuZYMOoa/zqRJjdc4NGTQZuLIGhVmb5preE6PrMph7sswKea+wsptJDfVmPnPU4r4sVcog46gis8UhYrZ7u2W3ARrDLYK6jiFPJI7a/Ldv+jlXibGqt6VKvdmiLapoWTBVbZT6OJYI2sRQpw1YmzLIwC1rleoETllQJxaETdEANlf33ozzUu8cq5pKMtb0827zOVeps+r+28Zr12iLlqja3/T5ISfYgJGM9SnkWcGyfU3kLsg0gGfU/QjbtkWpd45VTSUZK2pedj6zorhp8zlXjbEqsdv3oEldtEUF6jQ0fSZICTZgJGN9KQu0OsFRtq9U/wMERgOYafoRtsnBU+8cq5pKMlbUvMx/jvNx0/ZzDnWel2DfuQttUYFUBwYSbMBIxvpSZRqyao9adfw/taHhwGgAM6E/wiqhM+Cw2WUqydii5qWoP5r9DFP5nFM/Z6QtKlAl6QmVgNXZTyqBPYNkrC8hgzJkVp/gGUJVNICZkB9hnamllDvHqqaSjLmf3bwsiplFTVJfn3NZ85fyagraogIhByGWqbufBBswkrGQ+lrH1SSrL9pfgmcIVdEAZurOcIdaI5Ry51jVlJKxeXWuZCz7nNsmUgn2jbXQFhUIuTxnmZB9YiQkY6H0uY4r5FnAgFtBGsBMyNGsAQ+UNjLlZKyvQYkq9Qz4nNDdaYuWCnnhWpE+p0M7QjIWSt/ruELOjycepEVoAM8Itc5r6J1iXVNOxtzDjHqVxUyVmBr6SQBtUUOh+qa+pkM7RDIWSsrruELeNiMhU2wAu74CsmwQdQBhUcvUk7EyIUZTQ978NVVTbIuCCDVr09d0aIdIxkKJvY4r1BqzAZxB7JhaA9j2o2mzRmhAYVELydhyIUZTqw70Dzm+ptYWBRWqH2w7HRr5bJNkLJSY67hCrjEbwBnEjqk1gG0/mjYdXqiwSG10berJWNnn0XY0tcr7VY8lZVNri3oRMoEawDQmyVhIsdZxhVxjNqDFG1NrAOt8NEUfd9MOL9R3WaY2+jHlZCzkwvoqV0sONdGqYmptUSsxEqgBTGOSjMUSch1XkwQq5BRpJFNrAOt0jGVtVN3OMcRUU4qhNeVkLOb04diSs6m1RY3FTKD6uKqzhejJmKQbJT0n6aikgwveP1/SA/n7j0naX7bPZAK3r3VcTYJyBLe6mA3escRRiJGlLkbkQ4yiJNDenaXrGPKU2qM5bUc7myZUA2piKhtjW9SJVBOoBM4UoyZjks6R9IKkA5LOk/SEpGvntvmUpLvzx7dKeqBsv0kEbp/ruOq2biGnSCPaCd6xxFGoEa2yNqppu9N2fVEC7d1Zuo4hT6U9WqDN55HC+sOUjK0t6kyf68Dq7CuBM4TYydj1kh6Zef45SZ+b2+YRSdfnj8+V9JokW7bfJAK373VcdRKoFIcoGphpAEcRR6E6qVgjVANYI3uWrmPIU2mPFoiVUI2k+dllbG1RbSmuA2tyIVzEQYg2ydgetfcOST+deX4if23hNu7+G0m/lLQ3QN3dOn68/PX1demll6TTp7OfkrR/v7RnT/Zzayt7bd++xfuafb3qvqrub9bWVvG+0jCKOKoSMlUcPiytrOx+bWUle12q//HPWhYKZfWur0ubm9LqqmQm7d0r/c7vSLfdlkRYjSKG6pj/PFZXs+fr6+W/WyVWi2KlTfwNwOTiSFtb0saGdOxYlvIcO5Y9X/QfuqyROHRIOnVq9/unTmWvzysL4Dr72tnfbB9a5T9CKppmcTtF0sck3TPz/DZJX5zb5ilJV8w8f0HSJQv2tSFpW9L2vn37Oslca4m5jivkGUOKwxk5nTkbHUUchVxO0cVVjSEvDEglrLqIIU+xPWph0WfaZoAjlc8+pLG1RbWkug5sYEOwYpqyIzHXcYWcIk14gYdGNjXQZ4686OMvC4mQoZBKWHUdQ55Ke9RQUZzddVe7a4oGsiy1srG1RbX0fT+wkP1gQmInY+dKelHSVTqz2PHdc9t8WrsXOz5Ytt9kAjfWOq5U9xXYTAM4mjiaDZm9e7PS9mLcqvW2/dqbOlIJq65jyFNqjxpYFmddDnAMLVkbY1uU3Dqw0DNEiYmajGX16yZJP8mHag/lr31e0s354wskfVXZZcA/kHSgbJ/JNn7LgjvkUG/MffVIuy8nH1UclbUjoZOZKh/zmEfGvKMY8pTbowqaxlmsKzVjGV1bFHKYPlR/1CSoBpTVR0/GuihJNn59ruOKeUVKj9oEb5USM47K2p067VKV9ijE194sUlR3KmHVdQx5qu1RRU2TqlhXasYyurYoxXVgqQynd4RkrC99r+OKsa+eja4BnFHW7lTt7KpuV7XtrRMKVfL82GFFMrZcm6RqWSK+7HMfYp87urYoxXVgQ8zSayAZ60uqC25CfpVSz0bXAM4IlSPXSbJCJ05DaDtJxsqF/O9eJbkbQtzMG11blOI6sFSG0ztCMtaXkPPdMfeV0H+I0TWAM0L9M9fJtUNPKQ5hhINkrJ4qiVnbJagJNTGVja4tSnUdWCIDAV0gGetLyCtBYu4rodPW0TWAc0K0OyE+rqb7SChUCpGMnVEWb1WailAXngytzx1lW5TKOrChBUNDJGN9qhNUoebZQ+8roeGOUTaAgYUYZWj6kQ9hhINkLBNq+jDkhSdDMrm2qK91YENoRAIhGUtV3R4w5FlM23mGnkyuAWyo7Yllm/VrqZ/UkoxlqnzGVZqRUBeeDM3k2qK+1oEl1N90jWQsplBJT8iArrKSO5HWdHINYEfaTk8lFBK1kYxlqiRaIUbG3Ksn6Kkn8rMm2Rb1sQ4soZmYrpGMxRKyhws51BuyNe3YJBvAwKqGxkAGS2sjGcuEWlgfKjEfWoJPWzQjZAI1oP6oLZKxWEIGWchFkAM6E6EBbC9EIjWgkDnL1JKxtlfMtr2asuo2Q0vwJ9EWhVoHVrfOsQ7LzyEZiyXk/b2aBH/IfUUyiQYwkKKPO0QiNaCQOcuUkrEq/VofAwxV+s+iuNyJq9QGQUbfFoVcB9ak7jEOy88hGYulybx6ire6iGj0DWAgyz7SEG3ZgELmLFNKxlLpt9qsPZtP0lKJs9G3RSHXgYU05GH5OSRjsdTtwULeniLkviIafQMYyLKPO+QanwGEzFmmlIyFvMdXm8+76fegFv1eCoMgo2+L+k56YkyJRkYyFlOdFi32VyAlaPQNYE1NpyKHmkiFMKVkLJVF+lX7z/m4XPQ7qTRbo2+L+kx6Yk6JRkQylpKQc+MTmGcffQNYQ1dTkWNP1KaUjFXpt9pMIVZtPpr2nyk3W6Nvi/pMelKdEu0YyVgqQl41MpErUEbfANbQxVTkSMJkqSklY+7l/VaIG7uGOI6i30k1HifRFqWyDmwkydc8krFUhLzVRch9JWwSDWBFTduvCQygLjW1ZKxMHyNjy5Q1S6k2W7RFAfWxyDVBJGOpYE1YbTSAZzS9u0mIL3UeMpKxzOx0dtkVi131h0PuZ2mLAur68u9EtYmhPUI4+/bVe31rS9q/X9qzJ/u5tdV8Xxi8w4ellZXdr62sZK8XOXRIOnVq92unTmWvS4TRVGxtSRsb0rFj2XN3ySx7vLoqbW5K6+tntl9fz15bXc22W7RNE2XxiIlYFmDHjy/+ndnXl/WNY9U0i+u6DPIsgjVhtYmz0V3qTuFM9UudZ3UdQz6AOEplsGHII7G0RT2pclumgTZabWKIkbGQ6pxulp1CdnXqiqStr0svvSSdPp39LPu4y0a+CKNpqDLY0AdGYlGqbApgosOrJGOhzfem0uLh1iqtZ92eGZNTZWqTMBq/VJKgJlPtmJiyM8RUzix61ioZM7OLzeybZvZ8/vOigu1+a2Y/ysuRNnUOyuxCDvfs58ZG9noqrWcCiKPmFrVrt9+enUROabnF1GMolSRo6COxU4+jTs2uAzt0KAvORWeIU+0bm85vZtOj+gdJB/PHByX9fcF2v6q771HMr0/08t46JG0TR+FMMay6jiEfSBylesuIoaAt6tBE7sivWLe2kPScpMvzx5dLeq5gu2kG7kRvfFdH3gASR4GkspC7T13HkE8wjqaItqhDE7kjf5tkrO2ascvc/ZX88auSLivY7gIz2zazR83swy3rHI4qq6tZzCMRR8FMdLmFRAwhDOKoC3Ubpgn2jeeWbWBm35L0tgVv7bq0wd3dzLxgN6vu/rKZHZD0HTN70t1fWFDXhqQNSdo3hvnhw4ezNWKzV4ZMdDXrDTfcoFdffXXRW2+dfUIctbNv35l7Tc2/PnR9xpA07TgaM9qiCMbcMIXSdEjNa0xTzv3OvZI+VrbdaIZ0Bzrc2hfVmBrwKcdRRQNebtFY1zHkE4yjKaIt6tBEGiZFnKY8Iun2/PHtkr4+v4GZXWRm5+ePL5H0fknPtKx3OCY43NoAcRTI0K9ma4EYQgjEURcm3DBVVTpNWeILkh40s09IOibpLyXJzNYkfdLd75T0Lkn/r5mdVnYrjS+4O4GLWcRRQOvrk2zjiCGEQBx1ZaINU1WtkjF3f13SBxe8vi3pzvzxf0r6/Tb1YNyII7RFDCEE4gixcAd+AACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIiIZAwAACAikjEAAICISMYAAAAiIhkDAACIiGQMAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIiIZAwAACCiVsmYmf2FmT1tZqfNbG3Jdjea2XNmdtTMDrapE+NDHKEtYgghEEeIpe3I2FOS/lzS94s2MLNzJH1J0ockXSvp42Z2bct6MS7EEdoihhACcYQozm3zy+7+rCSZ2bLN3ivpqLu/mG/7FUm3SHqmTd0YD+IIbRFDCIE4Qix9rBl7h6Sfzjw/kb8G1EEcoS1iCCEQRwiudGTMzL4l6W0L3jrk7l8PeTBmtiFpI3/6hpk9FXL/NVwi6TXqDer3JL1lwevvDl1RInEU67OMWXfX9fYWQ9Lk42jM8UtbNP66Y9X7zqa/WJqMufsNTXeee1nSlTPPr8hfW1TXpqRNSTKzbXcvXEDZpVh1T63enborbjqoOIr9bzqlv7mLGJKmHUdTjd+Km9IWJV73AGLoLH1MU/5Q0tVmdpWZnSfpVklHeqgX40IcoS1iCCEQRwiu7a0tPmJmJyRdL+nfzeyR/PW3m9nDkuTuv5H0GUmPSHpW0oPu/nS7w8aYEEdoixhCCMQRYml7NeXXJH1twev/JemmmecPS3q45u432xxbS7Hqnlq9krQ50jgifnust+MYkib4bxqp3ph10xaNp+7B1WvuHvJAAAAAUANfhwQAABBRMslYzK+hMLOLzeybZvZ8/vOigu1+a2Y/ykvjBZtlf4OZnW9mD+TvP2Zm+5vWVbPeO8zs5MzfeGegev/ZzH5edFm3Zf4xP64fm9l7WtQVJY6mEkMV6yaOmtc7iTgihnZtN+gYyvdFHO1+v34cuXsSRdK7lN2j43uS1gq2OUfSC5IOSDpP0hOSrg1Q9z9IOpg/Pijp7wu2+1WAukr/BkmfknR3/vhWSQ/0VO8dkr7YwWf7J5LeI+mpgvdvkvQfkkzS+yQ9NrQ4mkIMEUfEEW0RMUQcdRNHyYyMufuz7v5cyWZvfg2Fu/9a0s7XULR1i6T78sf3SfpwgH0WqfI3zB7PQ5I+aLb8+zkC1dsJd/++pF8s2eQWSf/imUclvdXMLm9YV6w4mkIMVa27E8RRcLRFZyOG6iOOzlY7jpJJxirq6msoLnP3V/LHr0q6rGC7C8xs28weNbOmAV7lb3hzG88uo/6lpL0N66tTryR9NB9WfcjMrlzwfhf6/nqRLuqbQgxVrVsijpqaQhwRQ93W12cMScTRIrU/11a3tqjLevxqpTp1zz5xdzezoktMV939ZTM7IOk7Zvaku78Q+lgj+oakL7v7G2b2N8rOZP408jGdJVYcEUOVEUcN6519MvE4IoYa1jv7ZOIxJA0kjqSekzHv8auV6tRtZj8zs8vd/ZV8KPHnBft4Of/5opl9T9IfKZuzrqPK37CzzQkzO1fS70p6vWY9tet199k67lG29qAPdb+mJkocEUPV6iaOliOOiKGm9VWpt+cYkoijRWp/rkObpuzqayiOSLo9f3y7pLPOaMzsIjM7P398iaT3S3qmQV1V/obZ4/mYpO94viqwhdJ65+a0b1Z2d+k+HJH0V/kVKO+T9MuZYfYudBFHU4ihSnUTR61MIY6IoTOGHkMScbRI/TjywFcZNC2SPqJsXvUNST+T9Ej++tslPTyz3U2SfqIsgz8UqO69kr4t6XlJ35J0cf76mqR78sd/LOlJZVdsPCnpEy3qO+tvkPR5STfnjy+Q9FVJRyX9QNKBQH9nWb1/J+np/G/8rqRrAtX7ZUmvSPqf/DP+hKRPSvpk/r5J+lJ+XE+q4MqjlONoKjFEHBFHxBAxRByFjyPuwA8AABDR0KYpAQAARoVkDAAAICKSMQAAgIhIxgAAACIKkoxZj1++CgBAV+jPEEOokbF7Jd245P0PSbo6LxuS/ilQvQAAhHSv6M/QsyDJmPf45asAAHSF/gwx9LVmrO8vXwUAoAv0Zwiu1++mLGNmG8qGfXXhhRded80110Q+InTt8ccff83dL419HAAQGn3atLTpz/pKxip9aaa7b0ralKS1tTXf3t7u5+gQjZkdi30MAFBD5S+Bpk+bljb9WV/TlH1/+SoAAF2gP0NwQUbGzOzLkj4g6RIzOyHpbyW9RZLc/W5JDyv7Qs+jkk5J+usQ9QIAEBL9GWIIkoy5+8dL3ndJnw5RFwAAXaE/QwzcgR8AACAikjEAAICISMYAAAAiIhkDAACIiGQMAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIiIZAwAACAikjEAAICISMYAAAAiIhkDAACIiGQMAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAIKIgyZiZ3Whmz5nZUTM7uOD9O8zspJn9KC93hqgXAIDQ6NPQt3Pb7sDMzpH0JUl/JumEpB+a2RF3f2Zu0wfc/TNt6wMAoCv0aYghxMjYeyUddfcX3f3Xkr4i6ZYA+wUAoG/0aR3Z2pL275f27Ml+bm3FPqJ0hEjG3iHppzPPT+Svzfuomf3YzB4ysysD1AsAQGj0aR3Y2pI2NqRjxyT37OfGBgnZjr4W8H9D0n53/wNJ35R036KNzGzDzLbNbPvkyZM9HRoAALVMqk8LMaJ16JB06tTu106dyl5HmGTsZUmzZwVX5K+9yd1fd/c38qf3SLpu0Y7cfdPd19x97dJLLw1waAAA1EKfNiPUiNbx4/Ven5oQydgPJV1tZleZ2XmSbpV0ZHYDM7t85unNkp4NUC8AAKHRp80INaK1b1+916emdTLm7r+R9BlJjygLyAfd/Wkz+7yZ3Zxv9lkze9rMnpD0WUl3tK0XAIDQ6NN2CzWidfiwtLKy+7WVlex1BFoz5u4Pu/vvufv/6e6H89f+H3c/kj/+nLu/293/L3f/v939/wtR7xhxtQkAxEWfdkbVEa2yvmt9XdrclFZXJbPs5+Zm9joC3GcM4ezMze8MCe/MzUsELACgf4cP7+6XpLNHtKr2Xevr9GVF+DqkgNqOanG1CQAgJVVGtOi72iMZq2FZshXiihOuNgEApGZ9XXrpJen06ezn/OhWlb6LJTjLkYxVVJZshTgz4GoTAMDQlPVd3PC1HMlYRWXJ1rIzg6pnBFxtAgAYmrK+i2nMciRjFZUNwxadGVx8cfUzgjZXm1RJ+BgmBgCEVtZ3sQSnHMlYRWXDsEVnBlK9M4KyuflFqgwBM0wMAGir6KR+Wd/FEpxyJGMVlQ3DFp0Z/OIXi/c3f0bQZmSryhAww8QAgDaantSzBKcCd0+yXHfddZ6a++93X111N8t+3n9/+e+srrpnYbu7rK7u3u/Kyu73V1Z273/ZNmaL6zA78/tVtolB0rYnEG8UCoXSZUmxT6urSn9WpEn/OTRt+jNGxmrYGYb913/Nnt92W/naqypnBG1HtqoMATNMDABoo80tLJoswZkSkrGa6g7TVlmUXyXAl21TJeFjmBgA0Aa3sOgOyVhNTdZelZ0RtB3ZqpLw8b1gAIA2uIVFd0jGauriEt0QI1tVhoAZJgYANNX2FhbcXqkYyVhNTddeLQtCRrYAAF0JmQQ1vYUFU5glmq7877qkeuVJlSsfQ/zOVIirKSkUygRKrD6tz/5nWV1trsQcijb9GSNjNTUZoWIeHQAQQ5/9z7L+se4SnyajeUOeBiUZa2B+mFZaHgCh1pkNOdAAAP2r0//U6WPq3sKizhKfJlOag58GbTqk1nVJdZpyXpUh4LLh2So3w2s61Jz6jfbENCWFQplAidWnVZ0erNPHdL1cp8mUZgrToG36s+gBWlSGkoy1vcN+1QBtEmhDWKtGMkahUKZQUl8zVqePaZr4VB0caPKNMSl8y0yb/oxpypaqDAEvm0evOp/fZKqTtWoAMG1V1znX6WOaLr2pusSnyV0Lhv4tM0GSMTO70cyeM7OjZnZwwfvnm9kD+fuPmdn+EPWmoGoAFM2jVw3qJoHWxT3RAGDsxtanVbnHZJ0+psqd+MvWni1b41V2X81F+x/8t8w0HVLbKZLOkfSCpAOSzpP0hKRr57b5lKS788e3SnqgbL9DmaZsOxXYxXx+3X3HJKYpKRRKQmWqfVqoNWOhpkWLpjTL6o65RrpNfxYicK+X9MjM889J+tzcNo9Iuj5/fK6k1yTZsv2mGriLPuw2AVD3P0BRcNYN2lSQjFEolJTK1Pq0WfN9yV13FfdtRf1O1UGApmu8Uh5kiJ2MfUzSPTPPb5P0xbltnpJ0xczzFyTD35VHAAAgAElEQVRdsmy/KQZuV8lNl8lc7DOFMiRjFAolpTKlPm2Zpv1d1SSraVKVwkL9Im36s6QW8JvZhpltm9n2yZMnYx/OWbpaEN/mOyPLjonvowSAOFLs06reS6xpf1d17VnTNV5DX6hfJEQy9rKkK2eeX5G/tnAbMztX0u9Ken1+R+6+6e5r7r526aWXBji0sKouiC8L9pA3b2WRPgAENdo+rc6NUZv2LVWTrPV16fbbpXPOyZ6fc072fH7AYL6/vOmm4v0P+sboTYfUdoqy+fIXJV2lM4sd3z23zae1e7Hjg2X7TXFIt+09xaq8P7ufKtOLKc+fVyGmKSkUSkJlzH1ayHuJLeujQt3IvGibRWvZUlgj3aY/CxW8N0n6ibJ580P5a5+XdHP++AJJX5V0VNIPJB0o22cKgTsvxN32QyR0dY8pZSRjFAoltTLWPq3OeqsQV0wuS8qq9IV93Ig2pOjJWBclhcBdpCzjLwv2Kv8Z6gZV2zOUmEjGKBTKFEoKfVqoviXEoEKVvrBO8pjCwn6SsYSEGBkLFVRDGDUjGaNQKFMoKfRpVacGy07gQwwqMDK2uyR1NeUYlC1erLK4se632xctWOTrkAAAO8q+GqnqAv8qfVTZBQBV+sI6V1xO/g78XZUUziIWqbowcdk2Vd6vOh/fdhg4NjEyRqFQJlBS7dNmhfxGmKpTmW3706bbdqFNfxY9QItKioHb57RflaAKMQwcG8kYhUKZQonZp1VNUuou8A8xqDAmbfozpilr6HPar8rNWkMMAwMAxqvOvcXqLJGp0kf9zu+cebx37+4pUexGMlZD6Bustr1BXdl/nLL1AQCAcasziBDqBH4nAXx95ja4//3f9fYxNSRjNVQ5a6iaYNU5WylS5T/O/NmLNOA7FAMAaqkziBDqBJ6LxxpoOr/ZdRnimrE6c+Sh1nPVXdyY2hy+WDNGoVAmUGL1aTHWDg/h4rEutOnPGBmroeysoc7ZQJWzlSqjbHW+CJyzFQCYli7XDhf1UWP9Mu8ukYzVtCz5qTMcXBasIaYxqxzHstcBAMPW1drhZX3UsgRw0F/m3SGSsYDqnA2Una1UGcWqG9ScrQDA9NSZQalqWR9VlABK4QcZxoJkLKA6w8FlZytlo1hNRs641QUAIISyPmpRAshSmWIkYwHVHQ5edrZSNorVJKi51QUAIIQmMy0slSlGMtbAsunB9fVspGnfvizADh1qNgRbNorVNKi7GK4GAExLk5kWlsoUIxmrqWx6sM70YVlSt2wUi6AGAMRSZaZlvo+76SaWyhRqek+MrkuK9xlzD/d9kG3v+ZXiPcOaEPcZo1AoEyip9mldKeqj7ror7pd5d6lNf2bZ76dnbW3Nt7e3Yx/GWfbsycJqnlk29Vf2/o79+7NRs3mrq2fulF9mayubBj1+PBsRO3x4eNOOZva4u6/FPg4A6FKqfVpXQvRxQ9OmPzs39MGM3b59iwNsZ3qw7P0dIRYyrq8PL/kCAIwfi/XrYc1YTWWLFqsuamTNFwBgrOjj6iEZq6ls0WLV20dwzy8AwBBVueE4fVw9JGMNlN0eosrtI7jnFwBgaKreMWBZH8dXIp2t1QJ+M7tY0gOS9kt6SdJfuvv/XrDdbyU9mT897u43l+17aosdp4oF/ABSQZ9Wru3C/J1kbvam5Ssr4xiMaNOftR0ZOyjp2+5+taRv588X+W93/8O8lAYtAAAR0KeVaLswn69EWqxtMnaLpPvyx/dJ+nDL/QEAEAt9Wom2C/O5ynKxtsnYZe7+Sv74VUmXFWx3gZltm9mjZkZwAwBSRJ9Wou3CfK6yXKz0PmNm9i1Jb1vw1q5BRXd3MytagLbq7i+b2QFJ3zGzJ939hQV1bUjakKR9U/9kAADB0ae1s7Ouq+kNxw8fXrxmbOpXWbZdwP+cpA+4+ytmdrmk77n7O0t+515J/+buDy3bbiyLHbEcC/gBpII+rR9j+PaYRWIu4D8i6fb88e2Svj6/gZldZGbn548vkfR+Sc+0rBcAgNDo03pQ5fZPU9M2GfuCpD8zs+cl3ZA/l5mtmdk9+TbvkrRtZk9I+q6kL7g7gQsASA19GqJo9d2U7v66pA8ueH1b0p354/+U9Ptt6gEAoGv0aYiFO/ADAABERDIGAAAQEckYAABARCRjAAAAEZGMAQAAREQyBgAAEBHJGAAAQEQkYwAAABGRjAEAAEREMgYAABARyRgAAEBEJGMAAAARkYwBAABERDIGAAAQEckYAABARCRjAAAAEZGMAQAAREQyBgAAEBHJGAAAQEQkYwAAABG1SsbM7C/M7GkzO21ma0u2u9HMnjOzo2Z2sE2dAAB0gT4NsbQdGXtK0p9L+n7RBmZ2jqQvSfqQpGslfdzMrm1ZLwAAodGnIYpz2/yyuz8rSWa2bLP3Sjrq7i/m235F0i2SnmlTNwAAIdGnIZY+1oy9Q9JPZ56fyF8DAGBo6NMQXOnImJl9S9LbFrx1yN2/HvJgzGxD0kb+9A0zeyrk/mu4RNJr1NuLd0aqF8AETbBPi9m+T61Pa9yflSZj7n5D053nXpZ05czzK/LXFtW1KWlTksxs290LF1B2KVbdU6t3p+4Y9QKYpqn1abHb9yn9zW36sz6mKX8o6Wozu8rMzpN0q6QjPdQLAEBo9GkIru2tLT5iZickXS/p383skfz1t5vZw5Lk7r+R9BlJj0h6VtKD7v50u8MGACAs+jTE0vZqyq9J+tqC1/9L0k0zzx+W9HDN3W+2ObaWYtU9tXpj1w0AbxppnzbF9n1w9Zq7hzwQAAAA1MDXIQEAAESUTDIW82sozOxiM/ummT2f/7yoYLvfmtmP8tJ4wWbZ32Bm55vZA/n7j5nZ/qZ11az3DjM7OfM33hmo3n82s58XXdZtmX/Mj+vHZvaeEPUCQCyx+rS++7N8X/Rpu9+v36e5exJF0ruU3aPje5LWCrY5R9ILkg5IOk/SE5KuDVD3P0g6mD8+KOnvC7b7VYC6Sv8GSZ+SdHf++FZJD/RU7x2SvtjBZ/snkt4j6amC92+S9B+STNL7JD0WOx4pFAqlTYnVp/XZn1X9G+jTyvu0ZEbG3P1Zd3+uZLM3v4bC3X8taedrKNq6RdJ9+eP7JH04wD6LVPkbZo/nIUkftJLv5whUbyfc/fuSfrFkk1sk/YtnHpX0VjO7vI9jA4AuROzT+uzPJPq0RWr3ackkYxV19TUUl7n7K/njVyVdVrDdBWa2bWaPmlnTAK/yN7y5jWeXUf9S0t6G9dWpV5I+mg+rPmRmVy54vwt8vQiAKeqi7euzP5Po0xap/bm2urVFXdbj11DUqXv2ibu7mRVdYrrq7i+b2QFJ3zGzJ939hdDHGtE3JH3Z3d8ws79Rdibzp5GPCQCSFKtPoz+rbDB9Wq/JmPf4NRR16jazn5nZ5e7+Sj6U+POCfbyc/3zRzL4n6Y+UzVnXUeVv2NnmhJmdK+l3Jb1es57a9br7bB33KFt70IfGnysAxBKrT0uoP5Po0xap/bkObZqyq6+hOCLp9vzx7ZLOOqMxs4vM7Pz88SWS3i/pmQZ1VfkbZo/nY5K+4/mqwBZK652b075Z2d2l+3BE0l/lV6C8T9IvZ4bZAWCsuujT+uzPJPq0Rer3aaGvMmhxdcJHlM2rviHpZ5IeyV9/u6SH565S+ImyDP5QoLr3Svq2pOclfUvSxfnra5LuyR//saQnlV2x8aSkT7So76y/QdLnJd2cP75A0lclHZX0A0kHAv2dZfX+naSn87/xu5KuCVTvlyW9Iul/8s/4E5I+KemT+fsm6Uv5cT2pgiuPKBQKZSglVp/Wd39W9DfQp9Xr07gDPwAAQERDm6YEAAAYFZIxAACAiEjGAAAAIiIZAwAAiChIMtbJl2YCANAz+jPEEGpk7F5JNy55/0OSrs7LhqR/ClQvAAAh3Sv6M/QsSDLmfBE0AGAE6M8QQ19rxvgiaADAGNCfIbhev5uyjJltKBv21YUXXnjdNddcE/mI0LXHH3/8NXe/NPZxAEBo9GnT0qY/6ysZq/Slme6+KWlTktbW1nx7e7ufo0M0ZnYs9jEAQA2VvwSaPm1a2vRnfU1T8kXQAIAxoD9DcEFGxszsy5I+IOkSMzsh6W8lvUWS3P1uSQ8r+0LPo5JOSfrrEPUCABAS/RliCJKMufvHS953SZ8OURcAAF2hP0MM3IEfAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIiIZAwAACAikjEAAICISMYAAAAiIhkDAACIiGQMAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIgoSDJmZjea2XNmdtTMDi54/w4zO2lmP8rLnSHqBQAgNPo09O3ctjsws3MkfUnSn0k6IemHZnbE3Z+Z2/QBd/9M2/oAAOgKfRpiCDEy9l5JR939RXf/taSvSLolwH4BAOgbfRp6FyIZe4ekn848P5G/Nu+jZvZjM3vIzK4MUC8AAKHRp6F3fS3g/4ak/e7+B5K+Kem+RRuZ2YaZbZvZ9smTJ3s6NAAAaqFPQ1AhkrGXJc2eFVyRv/Ymd3/d3d/In94j6bpFO3L3TXdfc/e1Sy+9NMChAQBQC31aj7a2pP37pT17sp9bW7GPKI4QydgPJV1tZleZ2XmSbpV0ZHYDM7t85unNkp4NUC8AAKHRp/Vka0va2JCOHZPcs58bG9NMyFonY+7+G0mfkfSIsoB80N2fNrPPm9nN+WafNbOnzewJSZ+VdEfbelNBVg8A4zH1Pq1Phw5Jp07tfu3Uqez1qTF3j30MC62trfn29nbsw1hqJ6ufDaaVFWlzU1pfP7PNoUPS8ePSvn3S4cNn3oNkZo+7+1rs4wCALg2hT+vbnj3ZiNg8M+n06f6Pp602/Rl34G+hLKuvOgTL6BoAYGr27av3+piRjLVw/Pjy16sMwTJnDgBISV8DBIcPZ7NJs1ZWstenhmSshbKsvixZk5gzBwCko88BgvX1bFnP6mo2Nbm6unuZz5SQjLVQltVXGYKtkrABANCHLgYIlo20ra9LL72UrRF76aVpJmISyVgrZVl9lSFY5swBAH0pm4IMPUDAUpxqSMZaWpbVVxmCZc4cANCHKolR6AGCKhe6cQEbyVjnyoZgmTMHAPShyhRk6AGCZSNtjJqdQTKWAObMAQBdqzIFGXqAYNlIGxewnUEy1lCVodVQ24Q4FgDAtFWdggw5QLBspI0L2Ga4e5Lluuuu81Tdf7/7yop7NrCalZWV7PXQ24Q4lpRJ2vYE4o1CoVC6LCn0aYv6i/POc9+7193MfXW1vO+4//5su6rbL/ud1dXdx7JTVldb/JERtenPogdoUUkhcItUCaBQ27gvD/6hBzPJGIVCmUJJpU+b7U/27nV/y1u88sl86JP/oQ8mzGvTnzFN2UCVodVQ25QtcGSYFwBQ1ewU5P/6X9L//M/u95et2aqzxqvK8hkuYDuDZKyBKvPuobYpC37uUwYAaKLuyXzV1+tcJckFbBmSsQaqXPobapuy4Oc+ZQCAJuqezFd9PdZVkoO+mK3p/GbXJZX59SJVFjGG2KbKmrAmCypTIdaMUSiUCZQU+7SiNVt33bW4T6m6xsvMF/ZbZmf2E7rPSmH9WZv+LHqAFpUUAzeGPgIsZjJHMkahUKZQUu3T5tv/u+5a3udU6S+WDSJ01aelcDEbyVjP+k5e2tRX9ruxzyZIxigUyhRKyn3arBBJzbJ+pen+y/qystG4PpCM9ahq8pLC1GGVY419NkEyRqFQplBS7dPmLUtq6vRrRds2SZqG0Je5O8lYn6qu4QqRsLVN6Koca+yzCZIxCoUyhZJqnzavqN/YuzdMv9YkaQrZ73aJZKxHVZKXEIETIrBCHWuXSMYoFMoUSqp92ryivmfvXm/dr1XdZl7VQYPYM1IkYz0KNdpUtp++7s4f+2yCZIxCoUyhpNqnLbKoXwl5cl83aYo9aFBV9GRM0o2SnpN0VNLBBe+fL+mB/P3HJO0v22eqgRtq7rossKsEfqjRNa6mpFAolDNlSn1aVTGXvcQeNKgqajIm6RxJL0g6IOk8SU9IunZum09Jujt/fKukB8r2m3LghrhCMcTI2BjuQUYyRqFQUipT7NOqiL2IPvW+zN2jJ2PXS3pk5vnnJH1ubptHJF2fPz5X0muSbNl+xxC4bRK2KoEfe/F9CCRjFAolpUKfVizUrZKGkFg10aY/C/F1SO+Q9NOZ5yfy1xZu4+6/kfRLSXsD1J2ssu/bKvuC1CpfoBrieykH/fURABAefVqBnX7tX/81e37bbbv7jSr9VtXvrZxc39Q0i9spkj4m6Z6Z57dJ+uLcNk9JumLm+QuSLlmwrw1J25K29+3b11HuOh5t59FTmIcXI2MUCiWhQp+2XNt+YwgXljXVpj8LMTL2sqQrZ55fkb+2cBszO1fS70p6fX5H7r7p7mvuvnbppZcGOLSwqmTqdbL5tpl/lbOQZWJ9mSsAJGwyfVoTVfuNov7t+PHF+519fZJ9U9Msbqcomy9/UdJVOrPY8d1z23xauxc7Pli239Tm10PfPyWFzD+FNWdiZIxCoSRUptKnNdX2Sv8h3Iy8qTb9WajgvUnST5QN1R7KX/u8pJvzxxdI+qqyy4B/IOlA2T5TC9xQVzfW3bbtQse29yHrGskYhUJJrUyhT1umbb/R9ovCU+ibmoiejHVRUgvcKpl6nWw+xH3EyvRxl/+2SMYoFMoUSmp9WpE+rvQPdVVmakjGehBjZKzt6NkQ7kNGMkahUKZQUuvTioToN0KMbM3Xcddd5X3VkPuz6AFaVFIL3BhrxtqOng1h3p1kjEKhTKGk1qcVCdFvhL7fWFffeRkayVhPqgROncw8xNnFsm2GMO9OMkahUKZQUuzTFgnVb4Scigw9M9UVkrGRajt6lsKZQhmSMQqFMoUylD5tWb8RckCiTvIUes12V9r0ZyHuM4YFQtw9uO1d+NvehwwAMC1F/YZUfuf8qnfXl6rdb2xHlW+bCfGNNFE1zeK6LimfRaR0JcgQRr+WESNjFAplAiXlPq2KWLd3cmfNGIG7QKh7pIS86iP2FSRtkIxRKJQplJT6tCZ9RuipwrrJU+g1210gGetRiLsH173qcqiJVhUkYxQKZQollT6t6QhSF4vox9a/tenPWDNWU5V57rK56zrf7VV1/h0AgDJNv/fx8GFpZWX3aysr2et1tpm1vi699JJ0+nT2c8rrmUnGaqqySLAsIKsuXJzkl6UCADpTpf9ZdAFalQvC6lw0FuIit1FpOqTWdUllSHdeiJvZVR3KrXrT1yEP84ppSgqFMoGSSp9W1v80ncase4/N2Ivtu9CmP4seoEUllcDdMRtoe/dmpc2Xd1cJxK7+06SEZIxCoUyhpNKnlfUbTW6eWrcvSuEGrV0gGetYF0lP1StDQv+nSQ3JGIVCmUJJrU8r6n+a3Dy1bl+Uwg1au9CmP2PNWAVN1m6VzYdXWbhYNv9e56Z5AABIZ/c/0pn+ak9BVrDs5ql1+6LB36C1AyRjFdQNtJBXQS5L2ghoAEAb8/3Vb3979jbLroiU6vdFda+6nAKSsQrqBlqdkbQ2V5QQ0ACANhb1V5J0zjnVv0avyS0t+Kq+3c6NfQBDcPhwduYwG7DLAq3qSNrOGcnOfndG0KRqQbmzzaFD2b737cuOacoBDQCorqi/On06K1U06YvW1+mrZlm25iw9a2trvr29Hfsw3rS1VT3Q9u/PEqt5q6tn5ufrbDdmZva4u6/FPg4A6FJqfdoO+qFw2vRnTFNWVOdOwVWHbFmADwCIKfZyF27+miEZ60DV+XAW4AMAYoq5fouv/Duj1TSlmV0s6QFJ+yW9JOkv3f1/L9jut5KezJ8ed/eby/ad6pBuSPNrxqTsjGRKCxmZpgSQCvq0fo1tijTmNOVBSd9296slfTt/vsh/u/sf5qU0aKeCK0oAICn0aT1iqc4ZbZOxWyTdlz++T9KHW+5vcvjWegBIBn1aj1iqc0bbZOwyd38lf/yqpMsKtrvAzLbN7FEzI7gBACmiT+tR7IsHUlJ6nzEz+5akty14a9ctTN3dzaxoAdqqu79sZgckfcfMnnT3FxbUtSFpQ5L2TTE1BgB0ij4tHdwr84y2C/ifk/QBd3/FzC6X9D13f2fJ79wr6d/c/aFl27HYcRpYwA8gFfRpaCPmAv4jkm7PH98u6evzG5jZRWZ2fv74Eknvl/RMy3oBAAiNPg1RtE3GviDpz8zseUk35M9lZmtmdk++zbskbZvZE5K+K+kL7k7gAgBSQ5+GKFp9N6W7vy7pgwte35Z0Z/74PyX9fpt6AADoGn0aYuEO/AAAABGRjAEAAEREMgYAABARyRgAAEBEJGMAAAARkYwBAABERDIGAAAQEckYAABARCRjAAAAEZGMAQAAREQyBgAAEBHJGAAAQEQkYwAAABGRjAEAAEREMgYAABARyRgAAEBEJGMAAAARkYwBAABERDIGAAAQEckYAABARCRjAAAAEbVKxszsL8zsaTM7bWZrS7a70cyeM7OjZnawTZ0AAHSBPg2xtB0Ze0rSn0v6ftEGZnaOpC9J+pCkayV93MyubVkvAACh0achinPb/LK7PytJZrZss/dKOuruL+bbfkXSLZKeaVM3AAAh0achlj7WjL1D0k9nnp/IXwMAYGjo0xBc6ciYmX1L0tsWvHXI3b8e8mDMbEPSRv70DTN7KuT+a7hE0mvU24t3RqoXwARNsE+L2b5PrU9r3J+VJmPufkPTnedelnTlzPMr8tcW1bUpaVOSzGzb3QsXUHYpVt1Tq3en7hj1ApimqfVpsdv3Kf3NbfqzPqYpfyjpajO7yszOk3SrpCM91AsAQGj0aQiu7a0tPmJmJyRdL+nfzeyR/PW3m9nDkuTuv5H0GUmPSHpW0oPu/nS7wwYAICz6NMTS9mrKr0n62oLX/0vSTTPPH5b0cM3db7Y5tpZi1T21emPXDQBvGmmfNsX2fXD1mruHPBAAAADUwNchAQAARJRMMhbzayjM7GIz+6aZPZ//vKhgu9+a2Y/y0njBZtnfYGbnm9kD+fuPmdn+pnXVrPcOMzs58zfeGajefzaznxdd1m2Zf8yP68dm9p4Q9QJALLH6tL77s3xf9Gm736/fp7l7EkXSu5Tdo+N7ktYKtjlH0guSDkg6T9ITkq4NUPc/SDqYPz4o6e8LtvtVgLpK/wZJn5J0d/74VkkP9FTvHZK+2MFn+yeS3iPpqYL3b5L0H5JM0vskPRY7HikUCqVNidWn9dmfVf0b6NPK+7RkRsbc/Vl3f65ksze/hsLdfy1p52so2rpF0n354/skfTjAPotU+Rtmj+chSR+0ku/nCFRvJ9z9+5J+sWSTWyT9i2celfRWM7u8j2MDgC5E7NP67M8k+rRFavdpySRjFXX1NRSXufsr+eNXJV1WsN0FZrZtZo+aWdMAr/I3vLmNZ5dR/1LS3ob11alXkj6aD6s+ZGZXLni/C3y9CIAp6qLt67M/k+jTFqn9uba6tUVd1uPXUNSpe/aJu7uZFV1iuuruL5vZAUnfMbMn3f2F0Mca0Tckfdnd3zCzv1F2JvOnkY8JAJIUq0+jP6tsMH1ar8mY9/g1FHXqNrOfmdnl7v5KPpT484J9vJz/fNHMvifpj5TNWddR5W/Y2eaEmZ0r6XclvV6zntr1uvtsHfcoW3vQh8afKwDEEqtPS6g/k+jTFqn9uQ5tmrKrr6E4Iun2/PHtks46ozGzi8zs/PzxJZLeL+mZBnVV+Rtmj+djkr7j+arAFkrrnZvTvlnZ3aX7cETSX+VXoLxP0i9nhtkBYKy66NP67M8k+rRF6vdpoa8yaHF1wkeUzau+Ielnkh7JX3+7pIfnrlL4ibIM/lCguvdK+rak5yV9S9LF+etrku7JH/+xpCeVXbHxpKRPtKjvrL9B0ucl3Zw/vkDSVyUdlfQDSQcC/Z1l9f6dpKfzv/G7kq4JVO+XJb0i6X/yz/gTkj4p6ZP5+ybpS/lxPamCK48oFAplKCVWn9Z3f1b0N9Cn1evTuAM/AABAREObpgQAABgVkjEAAICISMYAAAAiIhkDAACIKEgy1smXZmJSiCGEQByhLWIIMYQaGbtX0o1L3v+QpKvzsiHpnwLVi/G4V8QQ2rtXxBHauVfEEHoWJBlzvggaLRFDCIE4QlvEEGLoa80YXwSNtoghhEAcoS1iCMH1+t2UZcxsQ9mwry688MLrrrnmmshHhK49/vjjr7n7pSH3SRxNSxcxJBFHU0NbhLbaxFBfyVilL810901Jm5K0trbm29vb/RwdojGzYxU3rfzFq8TRtNSIIYk4QgHaIrRVsy3apa9pSr4IGm0RQwiBOEJbxBCCCzIyZmZflvQBSZeY2QlJfyvpLZLk7ndLeljZF3oelXRK0l+HqBfjQQwhBOIIbRFDiCFIMubuHy953yV9OkRdGCdiCCEQR2iLGEIM3IEfAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIiIZAwAACAikjEAAICISMYAAAAiIhkDAACIiGQMAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIgoSDJmZjea2XNmdtTMDi54/w4zO2lmP8rLnSHqxbgQR2iLGEIIxBH61joZM7NzJH1J0ockXSvp42Z27YJNH3D3P8zLPW3rxbgQR2iLGEIIxFEAW1vS/v3Snj3Zz62t2EeUvBAjY++VdNTdX3T3X0v6iqRbAuwX00IcoS1iCCEQR21sbUkbG9KxY5J79nNjg4SsRIhk7B2Sfjrz/ET+2ryPmtmPzewhM7syQL0YF+IIbRFDCIE4auPQIenUqd2vnTqVvY5CfS3g/4ak/e7+B5K+Kem+RRuZ2YaZbZvZ9smTJ0QYu/AAABTHSURBVHs6NAwIcYS2KsWQRBxhKdqiIseP13sdksIkYy9Lmj0ruCJ/7U3u/rq7v5E/vUfSdYt25O6b7r7m7muXXnppgEPDgBBHaCtYDOXbTjqOJrzsh7aojX376r0OSWGSsR9KutrMrjKz8yTdKunI7AZmdvnM05slPRug3n5NuGXqyTTiCF0ihgKZ+LIf4qiqRf3i4cPSysru7VZWstdRqHUy5u6/kfQZSY8oC8gH3f1pM/u8md2cb/ZZM3vazJ6Q9FlJd7Stt1cTb5n6MIk46tjUzxeIoXCmvOyHOKqoqF+UpM1NaXVVMst+bm5mr0+5gSph7h77GBZaW1vz7e3t2IeR2b8/C7R5q6vSSy/1fTSjYmaPu/taV/tPKo46tNMuznagKytZG7i+Hu+4+tB1DEnTiaMde/Zk/es8M+n06f6Ppw+0RTXV6Rcn0kC1iSHuwF8FCxKRuCmPZCA8lv2gVJ1+kQaqFMlYFbRMCKjKdGLdKUfOFxASy35Qqk6/SANVimSsii5apqkv8JmoKssPmyxR5HwBO0I0Levri5f9jGhGCW3V6RdpoMq5e5Lluuuu86Tcf7/76qq7Wfbz/vurvVe0r5UV96yvzcrKSvnvjZCkbZ9QHK2u7v7Yd8rqar1t5k05pLqOIU8sjsqaoqnGQVtTa4uCWBSMRa9NIDDbxFBvyVXdEi1w+0ismvS2IzW1BtBs8UdvVm+bReqG7lhMKRkra276bFrGFm9Ta4s6sSxAxxYwC5CMhdJXYtW0tx2hqTWAXY2MTdmUkrGy2OgrkR/jQMfU2qJOTLzxahNDrBmb1eSKjyYLE7uYP2cN2iBUWWbB4mkUKWtumjQtTdYocnEcFmKhfmMkY7P6SqzKetu6iRU3pR2MKgujQy6eLgslcvhhKWtumiTyfZ2DYgJYqN9c0yG1rkuUId0+V04XzQtMbA2amBroTFkojWWqqesY8oTiqMpnVnfKscnU5oCbnEK0RTVNeKF+kTYx1GkD1qYMZs3Yzu+FutJyYmvQptoAVg2LKtsVbVMWSmPpUKeUjLmHXwe9LA5CnjOmbqptUSN1F+pPYPG+e7sY6rQBa1OSvJqyj1WuEztNnWIDWDUsqo6CFG1TFkoDzuF3mVoytiNUQl8UQ3fdVT6yOqb+dYptUWN1+py6/eCAA4tkrA99TR9O7CZTU2wAq37Eba+8ZGRsvHEUMqHf2W6+/xtLfFQ1xbaosTpncl0mbju/k0jyRjLWh76mD7tYg5ZIoC4yxQawali0vScZa8bGG0chE/oiYxk5rWpSbVHbfqFOYHWVuO38HQk1YiRjfehz+jDUVGligbrIpBrAXF8jY+7VpqgSztUrmWIyFjKhL8LI2EhjKMS0YZ19dJW41d13D0jG+tDF9GHXa9ASC9RFJtMAzuhrzdhUTDEZ62NkbGqxNZm2KNS0YdWF+l0lbu7JDd+SjPWhaWIV8nKkgQfqIpNpAOf0cTXlVEwxGQu9ZmxZPXUG4occh5Npi2Ks9+oicat7fD0gGetLzMTKffBDuItMpgFM3JA70ikmY+5hE/oQxzL0UbTJtEUprPcKkbiV7ScCkrHY+lrcP/DFjYtMpgEsEGJNV9tt2lxx1+RYQptqMlZXiM+v6P0qTVOq8bNjlG1Rquu9QiVuRX9jJCRjsfW1uL9JchXqYoCOjLIBrKhKGxNizVjZNlU70lTXr5GMZdom3G3isawJTDl+doyuLUp5vVeXV1hGRDIWW1+J1c7vhUiuYrd8udE1gDWUhU2VsAqxTZV2MdSxdIFkLEzC3SYe+4rlLo2uLUp5vVdXiVtkJGOxxU6smhxD7JYvN7oGsIYQd8gPsU2VUAh1LF0gGQuTcLeJx7LmJ+X42TG6tijl9V5dJW6RRU/GJN0o6TlJRyUdXPD++ZIeyN9/TNL+sn2m3vidJWZi5V4/YGO3fLnZ4J1aHKUyMlYl3FIe2eg6hjzxOHIPk3C3jcdlzVzK8bNjdG1RCtOGfSdukUVNxiSdI+kFSQcknSfpCUnXzm3zKUl3549vlfRA2X6TbfxSTKzcB3ul5U7wTi6OPJ01YzvbLAvrlNf8dB1DnngcuYdJuEPEY5GU42fH6NqiFKYNu0zcEhQ7Gbte0iMzzz8n6XNz2zwi6fr88bmSXpNky/abZOOXamLVpJ7YLV9upgGcThzNqJIElbVBobZpe6yh6qmr6xjygcRR24S7yjZtPt9U42fHKNui2NOGI12oXyR2MvYxSffMPL9N0hfntnlK0hUzz1+QdMmy/SbZ+KWaWLlXO62t8p+yZzMN4HTiCEF1HUM+kDhK4L/zoE2iLUp5vVciS2faaJOM7VFCzGzDzLbNbPvkyZOxD+dsx4/Xe12S9u2r97okHT4srazsfm1lJXu9yPq6tLkpra5KZtnPzc3s9a0taWNDOnYsC+9jx7LnkvTSS9Lp09nP9fXi/Q9I8nGEQRhaHK2vp//feWtL2r9f2rMn+7m11WyboUguhg4dkk6d2v3aqVPZ6/MBJBX3G4v6GunsD65OX9akrxyTplncTtGYhnTLdDFitez3Ql0QkPDwr8Y4NRAQ05Tluo4hH0EcVdX1NGXK68Ym0Ralst5r4Av1iyjyNOW5kl6UdJXOLHZ899w2n9buxY4Plu03icCdl0Ji1eQYEh7+nWkApxNHM8pCgwX85bqOIR9AHLm3T6S6XMDvnv4VlZNoi1JY7zWChfpFoiZjWf26SdJPlM2bH8pf+7ykm/PHF0j6qrLLgH8g6UDZPpMI3EViJ1ZNWqsBjIz51OLIyz/+UJ1XiCvtUu5Iu44hTzyO3MMkUmWfX5U4Wtb8pX6vsUm0RSms90q4P2orejLWRUkicOvoK7Fq0lolPPzbJnirlJTjqOzj56av1XQdQ554HLm3T6Tcu73pa9VjSGFkrKuSTAxVHTRIIXEbmDYxlNQC/kFbtjCySF8XBCxb3I9oyj7+Kh91iG2qhGGoY0E3yj7DEJ/xsverNH9V1nI3uXYJNVW90qOo35BYqN+Fpllc1yWZs4iqUr6FRcI0lbPRBUJMH4bYpkoYDmHNWJcl5ThyDzMy1maqs2rzl+pFIO7TbosqY6H+Um1iqNMGrE0ZXOD2mVgte31gQT7lBjDEwvoQ21QNm1Q7UpKxcIvvm14EMIZlQFNuiyrr+wrLgSEZS0EKidUAW8SpN4CptD+pHEcTJGOZtldTtq17YOeBZ5l6W1RJjCssB4RkLBWxE6sBLoykAVxuyElSX0jGyvUxqjn0WKUtmrPoA+UKy6VIxlLXV2I1wICeagNYtXMMNX04ZlNNxkJeFDeSgYlWptoWLVQUEHfdVfx6X4lbwkjGUpfK4v4ETbEBrPoxhVp0PfZEbYrJWOi7DgzwPC64KbZFhZYFxHyjsixB6yJxSxjJWOr6TKwG1vtOsQGsGg5VcviyNjPEBQKpm2IyFnrQoc754hhiZpEptkWFQk0xdpG4JYxkLHVdJFYjaRGn2ABWbeeqdLjL9lX2+wMcSF1oislY6OU4VZO7scTMIlNsiwrFWBtWNXFLONhIxoYgZGI1ohZxig1gyI5v2b5C3Hl/CKaYjIW+UC3k1PlQTbEtKhR6HnzHyBf1k4ylpI/EaoBBWmSKDWCdj7zK7QqK9lUWJiNZMzvJZKxusxHqasqxxMwiU2yLlqp6L7AUErdEkIylgltY1DbVBjDkLHPTO6qMJaefYjLmHmcVw1hiZpGptkWVLWtQQl7au2OAwUYylgpuYVEbDWC3yjrsMcx2TzUZK9L2c51CzCxCW1QiVL/TxYhbIkjGUsEtLGqbYgMYasqor2NJ3VSTsaLPrk2fOZUrcBeZYltUS1czMk1G3BINQpKxVHALi9qm1gCGXEyNzBSTsWUx0qbPbNqEjaApmlxbVFtXMzJ195twA0kylgpuYVHb1BrAkLcZQGaKydiyGGkTP3UTuYT7xdqm1hYt1ee0Yd2gS7iBbBNDe4Rw1telzU1pdVUyy35ubmavS9LWlrR/v7RnT/Zza+vM7730knT6dPZzdvuNDenYsSzcjh3Lnu/8Hgbn+PHy16tss0hReGF8lsXI4cPSysru11dWstfL7NtX7/VDh6RTp3a/dupU9joGqqjfkZb3b03VDbqmDWTqmmZxXZdBnUVUMfFbWBTRxM5GuxoZG9MIRV1dx5APMI6aDqjXjaMRXdg9ubaoUN/9Tt2gS7hfbBNDjIz1pckp5FjPACasyqhFk5ENRiimpSxGigbby5QN7s+rO6iBAei73ykKOmnxUH+bod+UNc3iui6DOYuoauK3sCiiCZ6NdnE15ZhGKOrqOoZ8wHHUxzGMZUR2im3RQin0O2WBlULwL9AmhloFl6SLJX1T0vP5z4sKtvutpB/l5UiVfQ8mcKti7mkhSdvEUXsptJ+xdB1DPqE4aiLRfrE22qJcCv3OQBu0NslY22nKg5K+7e5XS/p2/nyR/3b3P8zLzS3rHKayodVFq6/rzhkMF3HU0lhH7msghiJpOiWaKOIohX5nikt0mmZxWRKo5yRdnj++XNJzBdv9qu6+B3MWUUfT760ZMWVno8RRAGMZoair6xjyicXRVNEWJYSRsdouc/dX8sevSrqsYLsLzGzbzB41sw+3rHO4ik4hWX1NHAUwshGKuoghhEAcpWCCQ/3nlm1gZt+S9LYFb+3KFNzdzcwLdrPq7i+b2QFJ3zGzJ939hQV1bUjakKR9U7ocZwJDsjfccINeffXVRW+9dfYJcYQifcaQRByNFW3RAMwOVBw/nl2ee/jwqM8wLRtZa/jLZs9J+oC7v2Jml0v6nru/s+R37pX0b+7+0LLt1tbWfHt7u/GxDcr+/dmN9eatrmZDHCNmZo9L+j9EHKGhrmNIIo6mgLYIbZnZ4+6+1uR3205THpF0e/74dklfn9/AzC4ys/Pzx5dIer+kZ1rWOy4THJKdQxyhLWIIIRBHiKJtMvYFSX9mZs9LuiF/LjNbM7N78m3eJWnbzJ6Q9F1JX3B3AndWClevxEUcoS1iCCEQR4ii1TRllxjSnYY2w7pVEEfj13UMScTRFNAWoa2Y05QAAABogWQMAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIiIZAwAACAikjEAAICISMYAAAAiIhkDAACIiGQMAAAgIpIxAACAiEjGAAAAIiIZAwAAiIhkDAAAICKSMQAAgIhIxgAAACIiGQMAAIioVTJmZn9hZk+b2WkzW1uy3Y1m9pyZHTWzg23qxPgQR2iLGEIIxBFiaTsy9pSkP5f0/aINzOwcSV+S9CFJ10r6uJld27JejAtxhLaIIYRAHCGKc9v8srs/K0lmtmyz90o66u4v5tt+RdItkp5pUzfGgzhCW8QQQiCOEEsfa8beIemnM89P5K8BdRBHaIsYQgjEEYIrHRkzs29JetuCtw65+9dDHoyZbUjayJ++YWZPhdx/DZdIeo16g/o9SW9Z8Pq7Q1eUSBzF+ixj1t11vb3FkDT5OBpz/NIWjb/uWPW+s+kvliZj7n5D053nXpZ05czzK/LXFtW1KWlTksxs290LF1B2KVbdU6t3p+6Kmw4qjmL/m07pb+4ihqRpx9FU47fiprRFidc9gBg6Sx/TlD+UdLWZXWVm50m6VdKRHurFuBBHaIsYQgjEEYJre2uLj5jZCUnXS/p3M3skf/3tZvawJLn7byR9RtIjkp6V9KC7P93usDEmxBHaIoYQAnGEWNpeTfk1SV9b8Pp/Sbpp5vnDkh6uufvNNsfWUqy6p1avJG2ONI6I3x7r7TiGpAn+m0aqN2bdtEXjqXtw9Zq7hzwQAAAA1MDXIQEAAESUTDIW82sozOxiM/ummT2f/7yoYLvfmtmP8tJ4wWbZ32Bm55vZA/n7j5nZ/qZ11az3DjM7OfM33hmo3n82s58XXdZtmX/Mj+vHZvaeFnVFiaOpxFDFuomj5vVOIo6IoV3bDTqG8n0RR7vfrx9H7p5EkfQuZffo+J6ktYJtzpH0gqQDks6T9ISkawPU/Q+SDuaPD0r6+4LtfhWgrtK/QdKnJN2dP75V0gM91XuHpC928Nn+iaT3SHqq4P2bJP2HJJP0PkmPDS2OphBDxBFxRFtEDBFH3cRRMiNj7v6suz9XstmbX0Ph7r+WtPM1FG3dIum+/PF9kj4cYJ9FqvwNs8fzkKQPmi3/fo5A9XbC3b8v6RdLNrlF0r945lFJbzWzyxvWFSuOphBDVevuBHEUHG3R2Yih+oijs9WOo2SSsYq6+hqKy9z9lfzxq5IuK9juAjPbNrNHzaxpgFf5G97cxrPLqH8paW/D+urUK0kfzYdVHzKzKxe834W+v16ki/qmEENV65aIo6amEEfEULf19RlDEnG0SO3PtdWtLeqyHr9aqU7ds0/c3c2s6BLTVXd/2cwOSPqOmT3p7i+EPtaIviHpy+7+hpn9jbIzmT+NfExniRVHxFBlxFHDemefTDyOiKGG9c4+mXgM6f9v7+5Z4oiiOIw/BwKmk8Q0SWmVVkgR4jdIIQRSx8LG75AujZ/Azj6F3RaCEDV12iUpkmBlEAWLlGJxU8wVF992Znd2btZ9frDsMLvsmcv8i7Nz54UpyRF03IylDh+t1KR2RJxExPOU0nE+lHh6x2/8ye+HEfEVWKKas26izhguv3MUEY+AeeCsYZ3GdVNKgzW2qM496ELTx9QUyZEZqlfbHN3PHJmhUevVqdtxhsAc3abxfp22acpJPYaiB6zm5VXgxj+aiHgSEXN5+RmwDPwYoVadMQxuz3tgP+WzAscwtO61Oe0VqrtLd6EHfMhXoLwG/g4cZp+ESeRoFjJUq7Y5Gsss5MgMXZn2DIE5uk3zHKWWrzIY9QW8o5pXPQdOgN28/gWwM/C9t8BPqg7+Y0u1F4A94BfwBXia178CtvLyG6BPdcVGH1gbo96NMQCfgJW8/BjYBn4D34DFlsY5rO4G8D2P8QB42VLdz8AxcJH38RqwDqznzwPYzNvV544rj/7nHM1KhsyROTJDZsgctZ8j78AvSZJU0LRNU0qSJD0oNmOSJEkF2YxJkiQVZDMmSZJUkM2YJElSQTZjkiRJBdmMSZIkFWQzJkmSVNA/ApRsDHrCBUUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 12 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "clf = PCA(n_components=2)\n",
    "npos = [0, 0, 0, 0]\n",
    "npos = [smacof_mds(Cs[s], 2) for s in range(S)]\n",
    "\n",
    "npost01 = [0, 0]\n",
    "npost01 = [smacof_mds(Ct01[s], 2) for s in range(2)]\n",
    "npost01 = [clf.fit_transform(npost01[s]) for s in range(2)]\n",
    "\n",
    "npost02 = [0, 0]\n",
    "npost02 = [smacof_mds(Ct02[s], 2) for s in range(2)]\n",
    "npost02 = [clf.fit_transform(npost02[s]) for s in range(2)]\n",
    "\n",
    "npost13 = [0, 0]\n",
    "npost13 = [smacof_mds(Ct13[s], 2) for s in range(2)]\n",
    "npost13 = [clf.fit_transform(npost13[s]) for s in range(2)]\n",
    "\n",
    "npost23 = [0, 0]\n",
    "npost23 = [smacof_mds(Ct23[s], 2) for s in range(2)]\n",
    "npost23 = [clf.fit_transform(npost23[s]) for s in range(2)]\n",
    "\n",
    "\n",
    "fig = pl.figure(figsize=(10, 10))\n",
    "\n",
    "ax1 = pl.subplot2grid((4, 4), (0, 0))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax1.scatter(npos[0][:, 0], npos[0][:, 1], color='r')\n",
    "\n",
    "ax2 = pl.subplot2grid((4, 4), (0, 1))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax2.scatter(npost01[1][:, 0], npost01[1][:, 1], color='b')\n",
    "\n",
    "ax3 = pl.subplot2grid((4, 4), (0, 2))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax3.scatter(npost01[0][:, 0], npost01[0][:, 1], color='b')\n",
    "\n",
    "ax4 = pl.subplot2grid((4, 4), (0, 3))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax4.scatter(npos[1][:, 0], npos[1][:, 1], color='r')\n",
    "\n",
    "ax5 = pl.subplot2grid((4, 4), (1, 0))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax5.scatter(npost02[1][:, 0], npost02[1][:, 1], color='b')\n",
    "\n",
    "ax6 = pl.subplot2grid((4, 4), (1, 3))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax6.scatter(npost13[1][:, 0], npost13[1][:, 1], color='b')\n",
    "\n",
    "ax7 = pl.subplot2grid((4, 4), (2, 0))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax7.scatter(npost02[0][:, 0], npost02[0][:, 1], color='b')\n",
    "\n",
    "ax8 = pl.subplot2grid((4, 4), (2, 3))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax8.scatter(npost13[0][:, 0], npost13[0][:, 1], color='b')\n",
    "\n",
    "ax9 = pl.subplot2grid((4, 4), (3, 0))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax9.scatter(npos[2][:, 0], npos[2][:, 1], color='r')\n",
    "\n",
    "ax10 = pl.subplot2grid((4, 4), (3, 1))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax10.scatter(npost23[1][:, 0], npost23[1][:, 1], color='b')\n",
    "\n",
    "ax11 = pl.subplot2grid((4, 4), (3, 2))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax11.scatter(npost23[0][:, 0], npost23[0][:, 1], color='b')\n",
    "\n",
    "ax12 = pl.subplot2grid((4, 4), (3, 3))\n",
    "pl.xlim((-1, 1))\n",
    "pl.ylim((-1, 1))\n",
    "ax12.scatter(npos[3][:, 0], npos[3][:, 1], color='r')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
